September  2012 


ARBIC  ET  AL. 


1577 


Nonlinear  Cascades  of  Surface  Oceanic  Geostrophic  Kinetic  Energy 
in  the  Frequency  Domain* 

Brian  K.  Arbic 

Department  of  Earth  and  Environmental  Sciences,  University  of  Michigan,  Ann  Arbor,  Michigan 

Robert  B.  Scott 

Institute  for  Geophysics,  Jackson  School  of  Geosciences,  The  University  of  Texas  at  Austin,  Austin,  Texas, 
and  Department  de  Physique  et  LPO,  Universite  de  Bretagne  Occidental,  CNRS,  Brest,  France 

Glenn  R.  Flierl 

Department  of  Earth,  Atmospheric,  and  Planetary  Sciences,  Massachusetts  Institute  of  Technology,  Cambridge,  Massachusetts 

Andrew  J.  Morten 

Department  of  Physics,  University  of  Michigan,  Ann  Arbor,  Michigan 

James  G.  Richman  and  Jay  F.  Shriver 

Oceanography  Division,  Naval  Research  Laboratory,  Stennis  Space  Center,  Mississippi 

(Manuscript  received  24  August  2011,  in  final  form  31  March  2012) 

ABSTRACT 

Motivated  by  the  ubiquity  of  time  series  in  oceanic  data,  the  relative  lack  of  studies  of  geostrophic  turbulence  in 
the  frequency  domain,  and  the  interest  in  quantifying  the  contributions  of  intrinsic  nonlinearities  to  oceanic 
frequency  spectra,  this  paper  examines  the  spectra  and  spectral  fluxes  of  surface  oceanic  geostrophic  flows  in  the 
frequency  domain.  Spectra  and  spectral  fluxes  are  computed  from  idealized  two-layer  quasigeostrophic  (QG) 
turbulence  models  and  realistic  ocean  general  circulation  models,  as  well  as  from  gridded  satellite  altimeter  data. 

The  frequency  spectra  of  the  variance  of  streamfunction  (akin  to  sea  surface  height)  and  of  geostrophic  velocity 
are  qualitatively  similar  in  all  of  these,  with  substantial  variance  extending  out  to  low  frequencies.  The  spectral  flux 
U(oj)  of  kinetic  energy  in  the  frequency  co  domain  for  the  QG  model  documents  a  tendency  for  nonlinearity  to 
drive  energy  toward  longer  periods,  in  like  manner  to  the  inverse  cascade  toward  larger  length  scales  documented 
in  calculations  of  the  spectral  flux  fl(/c)  in  the  wavenumber  k  domain.  Computations  of  II(<x>)  in  the  realistic  model 
also  display  an  “inverse  temporal  cascade.”  In  satellite  altimeter  data,  some  regions  are  dominated  by  an  inverse 
temporal  cascade,  whereas  others  exhibit  a  forward  temporal  cascade.  However,  calculations  performed  with 
temporally  and/or  spatially  filtered  output  from  the  models  demonstrate  that  II (op  values  are  highly  susceptible  to 
the  smoothing  inherent  in  the  construction  of  gridded  altimeter  products.  Therefore,  at  present  it  is  difficult  to  say 
whether  the  forward  temporal  cascades  seen  in  some  regions  in  altimeter  data  represent  physics  that  is  missing  in 
the  models  studied  here  or  merely  sampling  artifacts. 
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1.  Introduction 

In  this  paper,  we  examine  the  spectra  and  spectral 
fluxes  of  surface  oceanic  geostrophic  flows  in  the  fre¬ 
quency  domain.  We  will  calculate  spectra  and  spectral 
fluxes  from  the  outputs  of  a  highly  simplified  (idealized) 
two-layer  quasigeostrophic  (QG)  turbulence  model 
(Arbic  and  Flierl  2004)  and  a  high-resolution  ocean 
general  circulation  model  run  in  a  realistic  near-global 
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domain  (called  here,  the  realistic  model;  Hurlburt  and 
Thompson  1980;  Shriver  et  al.  2007),  as  well  as  from 
gridded  satellite  altimeter  data  (LeTraon  et  al.  1998; 
Ducet  et  al.  2000).  We  are  motivated  in  part  by  the  fact 
that  oceanic  data  is  often  in  the  frequency  domain 
(taking  a  current  meter  as  an  example).  However,  as 
noted  by  Ferrari  and  Wunsch  (2010),  geostrophic  tur¬ 
bulence  remains  much  better  theoretically  developed  in 
the  wavenumber  domain  (taking  the  substantial  litera¬ 
ture  on  spectral  slopes  in  wavenumber  space  as  an  ex¬ 
ample;  e.g.  Vallis  2006).  As  far  as  we  know,  this  paper 
represents  the  first  comprehensive  look  at  idealized  QG 
turbulence  in  the  frequency  domain.  We  will  describe 
spectra  and  spectral  slopes  computed  from  idealized  QG 
models  and  compare  the  results  to  those  from  realistic 
ocean  models  and  oceanic  datasets.  Finally,  another 
primary  motivation  for  this  study  is  the  interest  in 
quantifying  the  contributions  of  intrinsic  nonlinearities 
in  oceanic  dynamics  to  oceanic  frequency  spectra.  As 
described  below,  we  will  quantify  nonlinear  contribu¬ 
tions  through  computations  of  spectral  fluxes.  Some 
important  related  issues  discussed  in  the  paper  include 
the  Taylor  (1938)  hypothesis  relating  frequency  and 
wavenumber  spectra,  the  sensitivity  of  QG  turbulence 
statistics  in  the  frequency  domain  to  bottom  friction,  and 
the  impact  of  the  spatial  and  temporal  filtering  inherent 
in  the  creation  of  gridded  satellite  altimeter  data  on  the 
computed  spectral  fluxes. 

We  will  adapt  some  of  the  language  and  results  of 
geostrophic  turbulence  theory,  which  is  generally  asso¬ 
ciated  with  the  wavenumber  domain,  to  the  frequency 
domain.  For  instance,  we  will  build  an  analog  in  the 
frequency  domain  for  the  diagnosis  of  the  inverse  cas¬ 
cade  to  larger  spatial  scales  commonly  done  in  the 
wavenumber  domain.  Recently,  the  spectral  fluxes  II(/:) 
of  kinetic  energy  in  wavenumber  k  space  for  surface 
ocean  geostrophic  flows  have  been  computed  from  sat¬ 
ellite  altimetry  data  of  sea  surface  height  (Scott  and 
Wang  2005). 1  Following  upon  this,  II(/c)  has  been  com¬ 
puted  from  both  idealized  QG  turbulence  models  (e.g., 
Scott  and  Arbic  2007)  and  realistic  ocean  general  cir¬ 
culation  models  (e.g.,  Schlosser  and  Eden  2007)  (see 
also  Tulloch  et  al.  2011).  All  three  types  of  computations 
demonstrate  that  there  is  an  inverse  cascade  of  kinetic 
energy,  driven  by  nonlinear  advection,  toward  larger 
spatial  scales  taking  place  at  length  scales  sufficiently 
large  [when  compared  to  the  first  baroclinic  (BC) 
mode  deformation  radius  Ld ].  Motivated  by  these 


1  In  this  paper,  k  denotes  isotropic  wavenumber  except  where 
noted;  in  a  few  places  in  section  2,  k  and  /  denote  wavenumbers  in 
the  x  and  y  directions,  respectively. 


studies,  here  we  investigate  spectral  fluxes  n(<u)  and 
spectral  transfers 

of  kinetic  energy  in  frequency  co  space,  as  diagnostics  of 
the  impact  of  nonlinearity  on  the  frequency  spectra  of 
surface  ocean  geostrophic  flows.  Our  analysis  of  non- 
linearities  in  the  frequency  domain  is  somewhat  similar 
to  the  diagnostics  employed  by  Sheng  and  Hayashi 
(1990a,b)  on  realistic  atmospheric  models  and  atmo¬ 
spheric  datasets.  Note  that  our  focus  is  on  geostrophic 
flows,  computed  from  derivatives  of  sea  surface  height 
(more  precisely  streamf unction)  via  the  geostrophic 
relations. 

In  analogy  to  the  spectral  fluxes  in  wavenumber  space, 
we  might  expect  that,  at  low  frequencies,  the  spectral 
fluxes  in  frequency  space  will  exhibit  cascades  toward 
yet  lower  frequencies  (longer  time  scales).  Indeed,  the 
Taylor  hypothesis  (Taylor  1938)  suggests  that  under 
certain  conditions  (e.g.,  turbulent  velocities  u'  that  are 
small  compared  to  a  mean  flow  velocity  U)  the  wave- 
number  and  frequency  domains  are  related  simply,  via 
the  formula  co  =  Uk.  In  this  case,  spectra  in  wavenumber 
and  frequency  space  would  have  the  same  slopes,  and 
phenomena  in  the  wavenumber  domain  (e.g.,  the  in¬ 
verse  cascade)  should  have  analogs  in  frequency  space. 
Further  motivation  arises  from  analysis  of  realistic  at¬ 
mospheric  models  and  atmospheric  datasets  in  the  fre¬ 
quency  domain,  which  demonstrate  an  energy  transfer 
from  high-frequency  motions  to  low-frequency  motions 
in  the  extratropics  (Sheng  and  Hayashi  1990a,b),  and 
from  experiments  on  acoustic  turbulence  in  superfluid 
helium,  in  which  an  inverse  energy  cascade  in  frequency 
space  takes  place  under  certain  conditions  (Ganshin 
et  al.  2008;  Efimov  et  al.  2009,  2010). 

Substantial  variance  often  extends  out  to  the  lowest 
frequencies  in  geophysical  records,  set  by  the  duration 
of  the  record  (e.g.,  Richman  et  al.  1977;  Wunsch  2009, 
2010,  and  references  therein),  and  there  is  great  interest 
in  examining  potential  causes  of  this  low-frequency 
variability.  An  ultimate  goal  of  the  work  begun  here  is  to 
understand  the  degree  to  which  oceanic  frequency 
spectra  in  the  geostrophic  regime  (as  measured  by,  e.g., 
current  meters  or  satellite  altimetry)  are  caused  by  in¬ 
trinsic  nonlinearities  in  the  oceanic  dynamics  (a  process 
emphasized  in  this  paper)  as  opposed  to,  for  instance, 
the  oceanic  response  to  atmospheric  wind  forcing  (e.g., 
Hasselmann  1976). 

Another  topic  we  will  explore  in  this  paper  is  the 
sensitivity  of  eddy  statistics  in  frequency  space  to  the 
strength  of  linear  bottom  Ekman  friction.  Several  recent 
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papers  have  demonstrated  that  the  length  scales  and 
other  statistics  in  stratified  geostrophic  turbulence  are 
highly  sensitive  to  the  strength  of  linear  bottom  Ekman 
friction  (Arbic  and  Flierl  2003,  2004;  Riviere  et  al.  2004; 
Thompson  and  Young  2006,  2007;  Arbic  et  al.  2007). 
Arbic  and  Scott  (2008)  showed  that  a  qualitatively 
similar  sensitivity  holds  when  the  bottom  friction  is 
quadratic.  Motivated  by  these  results,  here  we  briefly 
examine  the  sensitivity  of  statistics  computed  in  the 
frequency  domain  to  linear  bottom  Ekman  friction. 

We  will  also  examine  the  impact  of  filtering  a  geo¬ 
strophic  streamfunction  in  space  and  time  on  the  esti¬ 
mated  spectral  fluxes.  This  topic  is  important  because 
the  Archiving,  Validation,  and  Interpretation  of  Satel¬ 
lite  Oceanographic  data  (AVISO)  gridded  satellite  al¬ 
timeter  products  (Le  Traon  et  al.  1998;  Ducet  et  al.  2000) 
used  here  to  compute  spectral  kinetic  energy  fluxes  are 
constructed  by  interpolating  along-track  data  from  two 
altimeters,  in  space  and  time,  onto  a  V30  Mercator  grid. 
The  Ocean  Topography  Experiment  (TOPEX)/Poseidon 
and  Jason-1  altimeter  data  utilized  in  this  construction 
have  high  along-track  resolution  (about  6  km)  but  a 
large  distance  between  adjacent  tracks  (about  300  km; 
Chelton  et  al.  2001).  Because  of  the  large  along-track 
distance,  the  actual  feature  resolution  of  AVISO  data  is 
much  coarser  than  6  km  (Chelton  et  al.  2011).  The  in¬ 
terpolation  procedure  used  to  construct  the  gridded 
AVISO  data  acts  as  a  smoother  of  the  raw  along-track 
data,  in  both  space  and  time.  We  will  examine  the  impact 
of  such  smoothing  on  spectral  flux  estimates  by  com¬ 
paring  fluxes  computed  from  unfiltered  versus  filtered 
outputs  of  the  idealized  QG  models  and  realistic  ocean 
models. 

We  will  examine  the  consistency  of  slopes  of  stream- 
function  and  geostrophic  velocity  variance  in  the  fre¬ 
quency  domain  between  idealized  QG  models,  realistic 
ocean  models,  and  gridded  satellite  altimeter  data.  The 
slopes  of  spectra  in  frequency  and  wavenumber  space 
are  of  great  practical  and  theoretical  interest.  The  ad¬ 
vent  of  satellite  altimetry  has  made  it  possible  to  com¬ 
pute  spectral  slopes  on  global  scales  and  to  compare 
the  results  with  theoretical  predictions.  See,  for  instance, 
Le  Traon  et  al.  (2008)  and  Xu  and  Fu  (2011),  who  com¬ 
pute  wavenumber  spectra,  and  Hughes  and  Williams 
(2010)  and  Wunsch  (2009,  2010),  who  compute  fre¬ 
quency  spectra.  See  also  Hayashi  and  Golder  (1977) 
and  Hayashi  (1982),  who  discuss  the  space-time  spectra 
of  realistic  atmospheric  models  and  datasets.  As  noted, 
if  the  Taylor  hypothesis  holds,  the  slopes  in  frequency 
and  wavenumber  space  will  be  identical.  Therefore,  if 
the  Taylor  hypothesis  holds,  the  slopes  of  streamfunction 
and  geostrophic  velocity  variance  spectra  in  frequency 
space  will  differ  by  —2,  as  in  wavenumber  space  (because 


geostrophic  velocity  is  the  derivative  of  streamfunction 
and  a  spatial  derivative  in  physical  space  is  equivalent  to 
multiplication  by  k  in  wavenumber  space).  However, 
there  are  other  considerations,  such  as  Kolmogorov-like 
arguments  that  predict  a  -2  slope  for  an  energy  cascade 
in  frequency  space  versus  -5/3  in  wavenumber  space 
(Tennekes  and  Lumley  1972). 

The  paper  is  organized  as  follows:  In  section  2,  we 
present  results  from  the  QG  turbulence  model.  After 
discussing  the  model  limitations  and  governing  equa¬ 
tions  and  displaying  some  typical  flow  fields,  we  display 
frequency  spectra  and  the  sensitivity  of  the  QG  results  in 
frequency  space  to  the  strength  of  bottom  Ekman  fric¬ 
tion.  We  then  describe  the  theory  behind  spectral  fluxes 
of  kinetic  energy  in  frequency  space  and  demonstrate 
the  critical  role  that  the  spectral  kinetic  energy  flux  term 
plays  in  an  “integral  energy  budget.”  Next,  we  demon¬ 
strate  the  necessity  for  averaging  over  many  grid  points 
in  order  to  estimate  spectral  fluxes  in  frequency  space. 
Following  that,  we  demonstrate  the  dependence  of 
estimates  of  spectra  and  spectral  fluxes  in  frequency 
space  on  the  duration  and  sampling  frequency  of  the 
time  series  involved.  Section  2  concludes  with  a  dis¬ 
cussion  of  the  substantial  impact  that  filtering  in  space 
or  time  can  have  on  spectral  flux  estimates  made  from 
the  QG  model.  In  sections  3  and  4,  we  present  spectra 
and  spectral  fluxes  in  frequency  space,  computed  from 
the  realistic  high-resolution  numerical  ocean  general 
circulation  model  and  from  gridded  satellite  altimeter 
data,  respectively.  To  complement  the  filtering  exer¬ 
cise  done  with  the  QG  model  and  to  further  aid  in  in¬ 
terpreting  the  altimeter  results,  section  3  also  presents 
spectral  fluxes  computed  from  temporally  smoothed 
versions  of  the  realistic  model  output.  A  summary  and 
discussion  is  given  in  section  5. 

2.  Results  from  idealized  quasigeostrophic  model 

a.  Model  description,  governing  equations,  and 
example  flow  fields 

Our  two-layer  QG  turbulence  model  is  highly  simpli¬ 
fied.  There  is  no  planetary  beta  effect  or  bottom  topog¬ 
raphy,  and  there  are  only  two  vertical  modes  [barotropic 
(BT)  and  first  mode  baroclinic].  The  lack  of  planetary 
beta  excludes  Rossby  waves  (e.g.,  Pedlosky  1987;  Vallis 
2006),  which  are  an  important  component  of  surface 
oceanic  variability  (see,  e.g.,  Fig.  5  in  Wunsch  2010). 
Nevertheless,  we  will  see  that  the  frequency  spectra  and 
spectral  fluxes  in  the  idealized  QG  model  compare  well 
qualitatively  with  those  in  our  realistic  ocean  model. 

The  QG  model  is  forced  by  an  imposed  mean  flow, 
which  is  horizontally  homogeneous  within  each  layer  but 
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vertically  sheared  between  the  layers  and  is  therefore 
baroclinically  unstable  on  the  /  plane,  which  we  assume 
here  for  simplicity.  The  imposed  mean  flow  may  be 
thought  of  as  representing  the  time-mean  gyre  flow  in 
a  large  patch  of  the  ocean.  The  model  is  doubly  periodic 
and  has  256  grid  points  on  a  side  (thus  2562  grid  points  in 
all),  where  each  side  has  dimensional  length  207 rLd.  The 
model  is  initialized  with  a  randomly  generated  initial 
condition  and  achieves  statistical  equilibrium  when  the 
energy  extracted  by  eddies  from  the  mean  flow  is  bal¬ 
anced  by  energy  dissipation. 

As  in  our  earlier  papers  on  QG  turbulence,  the  gov¬ 
erning  equations  of  the  model  are  written  as  in  Flierl 
(1978).  The  Flierl  (1978)  equations,  which  can  be  written 
in  either  modal  or  layer  form,  bring  out  the  two  dy¬ 
namically  significant  properties  of  the  system-the  decay 
scale  as  embodied  in  a  single  internal  deformation  ra¬ 
dius  and  the  baroclinic  self-interaction  term.  In  layer 
form,  the  governing  equations  are 


_d<?l  , 
u,  — -  + 

1  dx 

^*1  and 

(4) 

dy  dx 

_ 3<72 

— -  + 

dq  2  dip2 

(5) 

2  dx 

~dy~dx ’ 

where  fluctuations  from  the  mean  are  written  without 
overbars,  whereas  imposed  time-mean  quantities  are 
designated  with  overbars.  The  zonal  velocity  is  u ,  the 
meridional  velocity  is  u,  and  the  imposed  mean  flow  is 
taken  to  be  zonal.  The  imposed  mean  PV  gradients  are 


dq1  (w1  —  u2)  dq2  S(u2  —  w1) 

~  (i  +  S)L2d'  ~  (i  +  s)l2 ' 


The  bottom  boundary  layer  thickness  dEkman  =  y/A/fQ 
determines  the  bottom  Ekman  drag  R2  by 


R 


2 


/p^Ekman 

2  H2  ’ 


(7) 


^  +  ri  +  =  ssd  and  (1) 

+  1*2  +  J(ip2,q2)  =  —R2y^*lf2  +  ssd,  (2) 

where  t  is  time;  potential  vorticity  q  is  the  prognostic 
variable;  ssd  indicates  small-scale  dissipation;  and  the 
subscripts  1  and  2  denote  the  upper  and  lower  layers, 
respectively.  The  perturbation  streamf unctions  i/q  and  ifj2 
satisfy 


<?1  =  V2^  + 


(■ft 2  ~  •h) 

(1  +  S)L2; 


q2  =  v+2  + 


SOi  -  «ft2) 

(1  +5)L2’ 


(3) 


where  8  =  HxlH^  the  ratio  of  upper-  to  lower-layer  depths, 
and  the  first  baroclinic  mode  deformation  radius  Ld  is 
defined  as  in  Flierl  (1978)  to  be  [g' H1H2lf^(H1  +  H2)]112, 
where  g'  is  reduced  gravity  and  /0  is  the  Coriolis  pa¬ 
rameter  (e.g.,  Vallis  2006). 2  The  upper-layer  stream- 
function  is  the  quantity  in  the  QG  model  that  is  most 
closely  analogous  to  sea  surface  height.  The  zonal 
(east-west)  and  meridional  (north-south)  coordinates 
are  %  andy,  respectively,  and /(A,  B)  =  dA/dx  dB/dy  - 
dA/dy  dB/dx.  The  forcing  terms  are 


2  Note  that  in  Vallis  (2006)  the  two-layer  QG  model  is  defined  to 
have  a  single  deformation  radius  as  is  the  case  here.  Some  treat¬ 
ments  of  the  QG  equations  differ  slightly.  For  instance,  in  Pedlosky 
(1987),  the  two-layer  QG  model  contains  two  distinct  internal  de¬ 
formation  radii. 


where  A  is  the  vertical  eddy  viscosity. 

The  ssd  used  here  is  an  exponential  cutoff  filter  in 
wavenumber  space  (Canuto  et  al.  1988;  LaCasce  1996), 


filter  =  exp 


/  _ _  \  M 

-a[Vk2  +  / 2  -  /c()J 


filter  =  1.0  when  Vk2  +  /2  </cn, 


when  Vk2  +  l2  >  k0 , 
(8) 


where  k  and  /  are  wavenumbers  in  the  %  and  y  directions, 
respectively,  and  /c0  is  the  cutoff  wavenumber.  We  take 
a  —  18.4,  M  —  4,  and  /c0  =  0.65 kN,  where  kN  corresponds 
to  the  Nyquist  scale.  The  filter  is  applied  to  qx  and  q2, 
the  Fourier  transforms  of  qi  and  q2 ,  at  every  time  step. 
As  shown  by  LaCasce  (1996),  the  wavenumber  filter 
does  not  distort  modons  (Stern  1975,  Larichev  and  Re¬ 
znik  1976a,b),  which  are  exact  solutions  of  the  inviscid 
nonlinear  QG  equations,  nearly  as  much  as  eddy  vis¬ 
cosities  or  hyperviscosities  do.  This  suggests  that  the 
wavenumber  filter  is  a  very  “clean”  choice  for  small-scale 
dissipation.  A  disadvantage  of  the  filter  is  that  it  is  not 
simple  to  diagnose  its  impact  on  the  energy  and  enstrophy 
budgets,  because  it  acts  on  q  instead  of  dq/dt.  We  infer  its 
impact  as  a  residual  of  the  other  terms  in  the  budget. 

BT  and  BC  modes  are  defined  as  in  Flierl  (1978), 


^BT 

9bt 


i  +  2 

1  +  8  ’ 
8qx  +  q2 
1  +  5 


,  _  -  <A2) 

'/'bc  i  +  a  > 

^  *AbT’ 


9bc 


 y/8(<h  ~  q2)  


^  ^BC  T  2^BC* 
^ d 


(9) 
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(a)¥l,FL=0.1  (b)q,,FL=0.1 


Fig.  1.  Snapshots  of  (left)  upper-layer  streamfunction  and  (right)  upper-layer  potential 
vorticity  qi  for  idealized  two-layer  QG  simulations  with  nondimensional  bottom  friction  values 
(top)-(bottom)  Fl  =  0.1,  0.4,  and  1.6. 


The  QG  equations  above  are  governed  by  two  non- 
dimensional  parameters.  In  the  runs  performed  here,  8  is 
set  to  0.2,  a  value  that  yields  a  strength  of  the  first  baro- 
clinic  mode  self-interaction  in  the  two-layer  model 
consistent  with  that  obtained  from  typical  midlatitude 
stratification  profiles  (Flierl  1978;  Fu  and  Flierl  1980). 
The  nondimensional  linear  bottom  Ekman  friction 
strength  FL  =  7?2Ly(Zq  -  uf)  is  set  to  0.4  in  the  “nom¬ 
inal”  run  presented  here,  which  we  will  focus  most  of  our 
efforts  on.  This  FL  value  yields  amplitudes,  horizontal 
scales,  and  vertical  structure  of  eddy  kinetic  energy  that 
agree  reasonably  well  with  midlatitude  observations  (Arbic 
and  Flierl  2004). 


Figure  1  displays  snapshots  of  i/q  and  qi  for  FL  values 
of  0.1  and  1.6  as  well  as  the  nominal  value  of  0.4.  The 
Fl  =  0.4  flow  field  consists  of  densely  packed  weak  vor¬ 
tices,  not  unlike  a  typical  snapshot  of  sea  surface  height 
in  the  midocean  (see,  e.g.,  comparisons  of  model  and 
observations  in  Fig.  10  of  Arbic  and  Flierl  2004).  In  the 
Fl  =  0.1  (weaker  friction)  experiment,  isolated  co¬ 
herent  vortices  (e.g.,  McWilliams  1984)  dominate  the 
i/q  and  especially  qi  fields. 

b.  Frequency  spectra  of  idealized  model 

As  is  standard,  the  model  output  is  windowed,  so  that 
the  beginning  and  ending  of  the  time  series  taper  toward 
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Fig.  2.  (a)  Wavenumber  k  spectrum  of  upper-layer  kinetic  energy  |Vi/q|2  in  two-layer  QG 
turbulence  model,  (b)  Frequency  co  spectrum  of  |Vi/q|2.  (c)  Wavenumber  spectrum  of  upper- 
layer  streamfunction  variance  i jj\.  (d)  Frequency  spectrum  of  ijj].  Cyan,  red,  and  blue  curves 
denote  experiments  with  nondimensional  bottom  friction  values  FL  =  0.1,  0.4,  and  1.6,  re¬ 
spectively.  Thick  solid  black  lines  have  been  included  to  show  slopes  of  (a)  -5/3  and  -3,  (b)  -2, 
and  (c)  -11/3  and  -5.  Thick  dashed  black  lines  in  (a)-(d)  are  least  squares  fits  to  the  red  curve 
over  the  range  of  wavenumbers  or  frequencies  covered  by  the  lines,  with  slopes  of  —3.7,  —1.8, 
-5.8,  and  -2.6,  respectively;  y  values  in  dashed  black  lines  are  multiplied  by  2  for  clarity.  In  this 
and  other  plots  of  QG  results,  U  in  the  x-axis  label  denotes  iq  -  zq. 


zero,  before  Fourier  transforms  are  computed  in  fre¬ 
quency  space.  We  utilize  nine  overlapping  Hanning 
windows,  each  having  width  equal  to  one-fifth  of  the 
total  time  interval.  Other  windows  were  tried  and  found 
to  yield  similar  results. 

The  wavenumber  and  frequency  spectra  of  the  upper- 
layer  streamfunction  variance  ifj\  and  of  the  upper-layer 
kinetic  energy  |Vi/q|2  are  displayed  in  Fig.  2.  The 


computations  in  Fig.  2  are  done  with  600  snapshots, 
where  each  snapshot  is  one  unit  of  nondimensional 
time  Ldl(U1  —  u2)  apart.  The  correlation  time  (as  mea¬ 
sured  by  the  zero  crossing  in  the  autocorrelation)  is 
about  16.5  snapshots  for  the  FL  =  0.4  run  and  a  similar 
number  for  the  other  runs  shown  in  Fig.  2. 

Results  in  Fig.  2  are  displayed  from  the  nominal  run 
with  Fl  =  0.4,  as  well  as  from  runs  with  FL  equal  to  0.1 
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and  1.6.  As  discussed  in  Arbic  and  Flierl  (2004)  and 
other  sources  (e.g.,  Thompson  and  Young  2006,  2007), 
as  the  friction  strength  decreases  {FL  decreases),  the 
inverse  cascade  reaches  out  to  larger  scales,  such  that  the 
wavenumber  spectrum  reddens  (cyan  curves  in  Figs.  2a, c). 
Figure  2a,  the  wavenumber  spectrum  plot  of  kinetic  en¬ 
ergy,  displays  the  canonical  slope  values  of  -5/3  and 
-3,  representing  the  classical  predictions  in  energy  and 
enstrophy  cascade  regimes  (Vallis  2006),  as  thick  solid 
black  lines.  Figure  2c,  the  wavenumber  spectrum  plot  of 
streamfunction  variance,  displays  the  classical  slope 
values  of  -n/3  and  -5,  obtained  from  the  -5/3  and  -3 
slopes  through  division  by  k2.  At  the  high-wavenumber 
end  but  before  the  sharp  dropoff  due  to  the  action  of 
the  wavenumber  filter,  the  slopes  of  the  streamfunction 
variance  and  kinetic  energy  in  the  FL  =  0.4  run  are  -5.8 
and  -3.7,  respectively,  as  determined  from  least  squares 
fits  over  the  range  indicated  by  the  thick  dashed  black 
lines.  These  slopes  differ  (to  within  “estimation  errors”) 
by  -2,  as  anticipated,  and  are  somewhat  steeper  than 
the  -5  and  -3  predictions  from  the  enstrophy  cascade. 

The  frequency  spectra  of  the  model  (Figs.  2b, d)  dis¬ 
play  a  steeper  slope  at  higher  frequencies  and  a  flatter 
slope  at  lower  frequencies.  This  basic  shape  is  also  found 
in  frequency  spectra  computed  from  current-meter  re¬ 
cords  (e.g.,  Richman  et  al.  1977;  Ferrari  and  Wunsch 
2010;  and  references  therein)  and  in  frequency  spectra 
computed  in  previous  studies  using  altimeter  data 
(e.g.,  Wunsch  2009,  2010,  and  references  therein).  In 
contrast  to  the  wavenumber  spectra,  which  turn  over  at 
low  wavenumbers,  the  frequency  spectra  continue  to 
be  energetic  (stay  flat)  at  the  lowest  frequencies.  In 
Fig.  2b,  the  thick  solid  line  is  at  a  slope  of  -2,  the 
predicted  slope  for  an  energy  cascade  in  the  frequency 
domain  (Tennekes  and  Lumley  1972).  In  Figs.  2b, d,  the 
thick  dashed  lines  represent  least  squares  fits  to  the 
frequency  spectra  of  velocity  and  streamfunction  var¬ 
iance,  respectively,  over  the  windows  covered  by  the 
lines,  from  the  FL  =  0.4  run.  The  slopes  over  these 
windows  are  —1.8  and  —2.6  for  velocity  and  stream- 
function  variance,  respectively.  The  velocity  variance 
slope  is  close  to  the  predicted  -2  slope,  as  are  slopes 
computed  from  current-meter-derived  frequency  spec¬ 
tra  (e.g.,  Richman  et  al.  1977;  Ferrari  and  Wunsch  2010; 
and  references  therein).  The  model  slopes  of  —1.8  and 
-2.6  are  not  equal  to  what  is  seen  at  the  high-wavenumber 
end  in  the  wavenumber  spectrum,  and  they  do  not  differ 
by  -2,  as  do  the  slopes  of  these  quantities  in  wave- 
number  space.  The  Taylor  (1938)  hypothesis  relating 
wavenumber  to  frequency  spectra  thus  fails  quantita¬ 
tively  in  both  respects.  At  the  low-frequency  end,  the 
slopes  of  upper-layer  streamfunction  and  velocity  variance 
are  -0.4  and  —0.1,  respectively,  as  determined  from 


least  squares  fits  over  a  decade  of  frequencies.  These 
slopes  also  do  not  differ  by  -2. 

c.  Sensitivity  of  time  scales  in  idealized  model  to 
Ekman  friction 

In  this  section,  we  compare  the  sensitivity  of  eddy  time 
scales  to  linear  bottom  Ekman  friction  with  the  previ¬ 
ously  documented  sensitivity  of  eddy  length  scales  to 
bottom  Ekman  friction.  As  is  common  (e.g.,  Vallis  2006), 
we  define  length  scales  C  as  reciprocals  of  the  first  mo¬ 
ment  (centroid)  of  the  wavenumber  spectrum  E(k).  For 
instance,  the  length  scale  £Bt  of  the  barotropic  mode  is 


where  EBT(k)  is  the  wavenumber  spectrum  of  baro¬ 
tropic  kinetic  energy.  Likewise,  we  define  model  time 
scales  T  as  reciprocals  of  the  centroid  of  the  frequency 
spectrum  E(co)  in  question.  For  instance, 


where  EBT(co)  is  the  frequency  spectrum  of  barotropic 
kinetic  energy.  For  these  quantities  to  be  well  defined, 
the  wavenumber  and  frequency  spectral  slopes  must  fall 
off  more  steeply  than  —2  at  the  high-wavenumber  and 
high-frequency  ends  and  less  steeply  than  - 1  at  the  low- 
wavenumber  and  low-frequency  ends.  The  nominal 
Fl  =  0.4  simulation  satisfies  these  constraints,  with  the 
exception  of  the  kinetic  energy  at  high  frequencies, 
which  falls  off  as  —1.8.  As  will  be  discussed  later,  a 
simulation  of  the  nominal  QG  solution  with  higher 
sampling  frequency  shows  slopes  steeper  than  -2  at 
the  high-frequency  end,  which  thus  satisfies  the  con¬ 
straints  for  coE  to  be  well  defined.  The  values  of  coE 
computed  from  the  two  cases  are  nearly  identical. 

The  length  and  time  scales  of  the  barotropic  mode, 
baroclinic  mode,  and  upper  layer,  in  QG  simulations 
with  different  bottom  friction  values,  are  displayed  in 
Fig.  3.  Length  scales  C  (dashed  curves  with  squares)  are 
normalized  by  Ld  while  time  scales  T  (solid  curves  with 
x’s)  are  normalized  by  Ldl(u1  —  u2).  In  the  weak  friction 
regime  (FL  of  order  one  or  less),  as  friction  further 
weakens  the  length  scales  of  the  barotropic  mode  and 
upper  layer  increase,  whereas  the  baroclinic  length  scale 
remains  relatively  constant.  In  the  weak  friction  regime, 
the  kinetic  energies  associated  with  each  mode  or  layer 
increase  more  rapidly  with  further  decreases  in  friction 
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(a)  Barotropic  T and  C 


(b)  Baroclinic  T and  C 


(c)  Upper  layer  T and  C 


Fig.  3.  Length  scales  C  and  time  scales  T  of  the  kinetic  energy  in 
the  (a)  barotropic  mode,  (b)  baroclinic  mode,  and  (c)  upper  layer,  in 
two-layer  QG  turbulence  experiments  damped  by  bottom  friction, 
as  a  function  of  nondimensional  bottom  friction  strength.  Length 
(time)  scales  are  computed  as  the  reciprocal  of  the  first  moment  of 
the  appropriate  wavenumber  (frequency)  spectrum.  Length  scales 
are  given  by  dashed  curves  with  squares  and  are  normalized  by  Ld , 
whereas  time  scales  are  given  by  solid  curves  with  x’s  and  are  nor¬ 
malized  by  Ldl(jl1  -  u2).  Dotted  curves  with  circles  are  the  predicted 
time  scales  obtained  by  dividing  the  length  scale  C  by  the  rms  ve¬ 
locity  scale  associated  with  the  layer  or  mode  in  question.  Dashed 
horizontal  lines  are  drawn  at  the  value  of  1. 

than  the  length  scales  do  (see,  e.g.,  Figs.  4-6  of  Arbic  and 
Flierl  2004).  Therefore,  if  we  crudely  estimate  the  time 
scale  by  T  =  £/Urms ,  where  t/rms  is  the  rms  velocity  of 
the  mode  or  layer  in  question,  then  we  might  expect 
a  decrease  in  T  with  decreasing  friction.  Figure  3  dem¬ 
onstrates  that  these  predicted  time  scales  (dotted  curves 
with  circles)  do  indeed  decrease  with  decreasing  friction 
strength.  The  actual  time  scales  computed  from  (11)  also 
decrease  with  decreasing  friction.  Thus,  the  trend  of  the 
changes  in  time  scale  with  changes  in  friction  is  suc¬ 
cessfully  predicted  by  the  above  argument.  However, 
the  actual  time  scales  are  larger  than  the  predicted  time 
scales.  The  actual  barotropic  time  scale  exceeds 
Ldl(uf  ~  u 2)  in  °f  the  experiments  analyzed  here.  In 
most  of  these  experiments,  the  time  scale  exceeds  this 


value  for  the  upper  layer  and  baroclinic  modes,  as  well  as 
for  the  barotropic  mode.  This  tendency  for  time  scales 
to  exceed  Ldl(u q  -  uf)  mimics  the  tendency  for  length 
scales  to  exceed  Ld ,  although  the  tendency  is  not  as  rigid 
in  the  case  of  the  time  scales.  Some  of  the  time  scales 
drop  below  the  dashed  lines  set  at  unity  in  Fig.  3,  whereas 
the  length  scales  never  drop  below  these  lines. 

d.  Theory  of  spectral  fluxes  in  frequency  space 
To  develop  the  theory  of  II(ca),  we  let 

k  l  oj 

<A 2{x,y,t)  =  1 1 1 (12) 

k  l  (o 

where  f  denotes  a  Fourier  transform  of  f  in  wavenumber- 
frequency  space  and  k  and  /  are  zonal  and  meridional 
wavenumbers,  respectively.  If  we  neglect  the  small-scale 
dissipation,  the  Fourier  transformed  upper-  and  lower- 
layer  governing  equations  are 


-icol  Kz  + 


(1  +  S)L 2, 


f/q(/c,  /,  co) 


+  i(o 


+ 


— —— — j  (k,l,to)  +  J(i//1,V2t//1)(k,  /,  co) 

(i  +  8)Ld 

JO/ThXk  l 0J)  +  L  ik,  /,&>)  =  0  (13) 


(1  +  8)L2 


and 


-i(o\K 


(1  +  8)0, 


8  W,  ,  .  .  8 

i/j2{kJ,(o)  +  Ko 


(1  +  $)L2d 


l/q(/c,/,<u) 


+  /(t^2,  V  ifj2)(k,l,(o)  +  ^  + 


+  r2(/c,/,(u)  —  i?2A^^2  (&,/,&>), 


(14) 


respectively,  where  K2  =  k2  +  l2  and  we  have  dropped 
the  summation  signs  for  simplicity. 

To  obtain  the  depth-averaged  energy  equation,  we  mul¬ 
tiply  (13)  by  -(5/1  +  ( k ,  /,  co ),  where  the  superscript 

asterisk  represents  a  complex  conjugate,  and  the  complex 
conjugate  of  (13)  by  —  (5/1  +  5)*/q (kj,  co);  add  the  re¬ 
sults;  and  divide  by  2.  Similarly,  we  multiply  (14)  by 
-(1/1  +  8)i/jf(kJ,co)  and  the  complex  conjugate  of  (14) 
by  -(1/1  +  8)iJj2(kJ,(o ),  add  the  results,  and  divide  by  2. 
We  then  add  everything  together  and  find  that  the  terms 
involving  ico  in  (13)  and  (14)  drop  out.  Utilizing  the  identity 
(a*b  +  ab*)/2  =  Re[a*5],  where  Re  denotes  the  real  part  of 
a  complex  number,  we  are  left  with 
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Re 


1+5 

Re 

Re 


+  Re 


1 


1  +  8 


if/2  (k,l,co)J(if/2,Wzil/2)(kJ,co) 


8 

§  -  —  1 

(1  +d)2L2^1  ^  (fc’/’wPW'i’^2)(M,w) 

+  Re 

1  +  8 


l/j2*  ( k ,  /,  CD)r2(/c,  /,  c o ) 


1 


1  +  5 


-R2K2ilj2*(kJ,(o)ilf2(kJ,(o)  =  0. 


(15) 


To  simplify  the  interpretation  of  (15),  we  utilize  the 
following  definitions:  The  spectral  transfers  of  upper- 
and  lower-layer  kinetic  energy  are 


TKB1(k,l,co)  =  Re  -  —  ^*(k,l,(D)J(ijjv\  j/q)(fc,/,w) 
and  (16) 

Tkb,2^  ^w)  =  Re 


=  Re 


(l+S)2L2d 


(*Ai  “  ik2)*(k,l,(x))J(ipv4f2)(k,l,a)) 


forcing  (/c’/’w)  =  Re 


+  Z,  o»)S^(*r,  /,  o») 


range  of  wavenumbers  yields  spectral  fluxes  in  wave- 
number  space, 


(17) 


respectively.  The  spectral  transfer  of  available  potential 
energy  (APE)  is 

^APE^’  ^?<w) 


nKE,lW 

^KE,2^)  =  | 

h 

nAPEW  = 

(K)=  [ 
J) 


k2+l2>K2  J 


k2+l2>K2 


k2+l2>K2 


forcing 


k2+l2>K2 


u 


friction  C^") 


k2+l2>K2 


TKE  x  (&,  /,  w)  dco  dk  dl ,  (22) 

2(^5  /,  w)  dco  dk  dl,  (23) 

Tape  (&,  /,  w)  dco  dk  dl ,  (24) 

^forcing (^»  U<*>)  dco  dk  dl,  and 
(25) 

^friction (^  ’  ^  w)  ^  dk  dk  (26) 


(18) 


Similarly,  the  spectral  fluxes  in  frequency  space  II(co)  are 
obtained  by  integrating  over  all  wavenumbers  and  over 
a  semi-infinite  range  of  frequencies, 


The  spectral  transfer  due  to  the  imposed  mean  flow 
forcing  is 


>  (19) 


and  the  spectral  transfer  due  to  bottom  Ekman  friction  is 

Th^JkJ^)  =  ~^R2K%*(k,l,<o)^2(k,l,<o). 

(20) 

With  these  definitions,  (15)  becomes 

TKE1(k,l,io)  +  TKE2(k,l,co )  +  l'APE(k,l,oj) 

+  7’fordng(/c’/’w)  +  7friction(/c’,’w)  =  °-  (21) 

Spectral  fluxes  are  defined  by  integrating  the  spectral 
transfers  over  some  part  of  (/c,  /,  co)  space.  For  instance, 
integration  over  all  values  of  co  and  over  a  semi-infinite 


iW«)  = 

nKE.2(-Q)  = 

nAPE(H)  = 
Hforcing(H)  = 

Hfriction(H)  = 


w>fl 


0>fl 


TKE 1  ( k ,  /,  co)  dco  dk  dl,  (27) 


TKE  2  (k,  /,  co)  dco  dk  dl,  (28) 


w>fl 


TApE(/c,  /,  co)  dco  dk  dl,  (29) 


w>fl 


0>fl 


^forcing  (^’  ^  <*>)  ^  dk  d ^  and 

(30) 

^frictioiA’  Z’ ")  dk  dL  (31) 


Following  the  terminology  of  Scott  and  Arbic  (2007),  we 
will  call  the  budget  obtained  by  summing  the  spectral 
fluxes  (22)-(26)  or  (27)-(31)  the  integral  energy  budget. 
Note  that,  because  we  are  using  the  depth-averaged  en¬ 
ergy  equation,  there  are  no  interlayer  spectral  fluxes  or 
transfers.  See  Sheng  and  Hayashi  (1990a,b)  and  refer¬ 
ences  therein  for  the  derivation  of  energy  equations  in  the 
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frequency  domain  for  realistic  atmospheric  models  and 
atmospheric  datasets. 

For  simplicity,  in  this  paper,  we  examine  either  spec¬ 
tral  fluxes  U(k)  in  wavenumber  space  or  II(m)  in  fre¬ 
quency  space.  In  future  work,  we  will  consider  fluxes  in 
which  partial  integrations  are  performed  over  both  fre¬ 
quency  and  wavenumber,  via 


n  (*,n)=  f 

Jy 


k2+l2>K2 


0)>fl 


T (k,  /,  co)  dco  dk  dl.  (32) 


In  the  next  section,  we  will  show  integral  potential  ens- 
trophy  budgets  as  well  as  integral  energy  budgets.  The 
potential  enstrophy  budget  is  obtained  through  multipli¬ 
cation  of  (13)  and  (14)  by  8q*/(l  +  8 )  and  qf/(  1  +  5), 
respectively,  followed  by  addition,  taking  the  real  parts,  and 
integration. 

e.  Energy  and  enstrophy  budget  of  idealized 
model  in  frequency  space 

The  integral  energy  budget  in  wavenumber  space  of 
the  nominal  FL  =  0.4  QG  simulation  was  examined  in 
detail  in  Scott  and  Arbic  (2007)  and  is  discussed  briefly 
here  (Fig.  4a).  The  forcing  term  (magenta  curve)  is  pos¬ 
itive,  indicating  an  energy  gain  over  all  wavenumbers. 
The  bottom  friction  term  (blue  curve)  is  a  sink  of  energy 
over  all  wavenumbers.  The  flux  nApE(k)  of  available 
potential  energy  (green  curve)  is  positive,  meaning  that 
potential  energy  cascades  toward  larger  wavenumbers 
(smaller  length  scales;  Salmon  1980).  The  black  and  red 
curves  represent  spectral  fluxes  nKE,i(&)  and  nKE,2(^)  of 
upper-  and  lower-layer  kinetic  energy,  respectively,  with 
the  latter  being  much  smaller  in  magnitude.  The  nega¬ 
tive  values  of  these  curves  indicate  an  inverse  cascade  of 
energy  toward  lower  wavenumbers  (larger  length  scales). 
The  cyan  curve  is  the  residual  of  the  terms  above.  It  is 
close  to  zero  over  all  wavenumbers,  indicating  that  the 
effects  of  the  wavenumber  filter  in  the  energy  budget  of 
these  runs  are  negligible,  consistent  with  our  results  in 
Arbic  and  Flierl  (2004). 

The  integral  energy  budget  in  frequency  space  is  dis¬ 
played  in  Fig.  4b.  The  terms  behave  similarly,  with  a 
positive  forcing  term  (energy  source),  negative  bottom 
friction  term  (energy  sink),  positive  potential  energy 
flux  (forward  cascade  to  small  time  scales),  negative 
kinetic  energy  flux,  and  small  residual.  The  negative  spec¬ 
tral  kinetic  energy  fluxes  IIke,^)  and  nKE?2(<y)  indicate 
that  an  “inverse  temporal  cascade”  exists  alongside  the 
well-known  inverse  cascade  in  space.  An  important  dif¬ 
ference  between  the  budgets  in  frequency  versus  wave- 
number  space  is  that  the  fluxes  in  the  frequency  budget 
are  “wider.”  The  forward  potential  energy  flux  and 


Fig.  4.  Integral  energy  budgets  of  nominal  two-layer  FL  =  0.4 
QG  simulation  in  (a)  wavenumber  k  space  and  (b)  frequency  co 
space.  Black  and  red  curves  denote  spectral  fluxes  of  upper-  and 
lower-layer  kinetic  energy  nKE,i  and  nKE,2>  respectively.  Green 
curves  denote  spectral  flux  of  APE  nAPE.  Forcing  is  given  by  the 
magenta  curve,  whereas  bottom  Ekman  friction  is  given  by  the  blue 
curve.  The  cyan  curve  denotes  the  residual  of  the  terms  above.  All 
terms  nondimensionalized  by  (zq  -  uf)3 1  Ld. 

inverse  kinetic  energy  flux  take  on  nonnegligible  values 
immediately  to  the  right  of  the  low-frequency  endpoints. 
In  contrast,  in  wavenumber  space  these  fluxes  are  close 
to  zero  over  wide  bands  at  the  low  end  of  the  spectrum. 
Nonlinearities  in  the  QG  turbulence  model  apparently 
drive  energy  to  the  lowest  frequencies  allowed  in  the 
simulation,  where  these  lowest  frequencies  are  set  by  the 
simulation  duration. 

The  integral  potential  enstrophy  budget  is  displayed 
in  Fig.  5,  in  both  the  wavenumber  and  the  frequency 
domains.  Once  again,  forcing  is  denoted  by  the  magenta 
curves  and  is  positive,  indicating  a  source  of  enstrophy. 
The  bottom  friction  terms  (blue  curves)  represent  a  sink, 
but  one  that  is  much  less  substantial  than  in  the  energy 
budgets.  The  flux  of  potential  enstrophy  by  nonlinear 
terms  (black  and  red  curves  for  the  upper  and  lower 
layers,  respectively)  is  positive,  indicating  a  forward  flux 
of  potential  enstrophy  down  to  small  spatial  scales  (where 
enstrophy  dissipation  takes  place)  and  small  time  scales. 
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Fig.  5.  Integral  potential  enstrophy  budgets  of  nominal  two- 
layer  Fl  =  0.4  QG  simulation  in  (a)  wavenumber  k  space  and  (b) 
frequency  co  space.  Black  and  red  curves  denote  spectral  fluxes  of 
upper-  and  lower-layer  potential  enstrophy,  respectively.  Forcing  is 
given  by  the  magenta  curve,  whereas  bottom  Ekman  friction  is 
given  by  the  blue  curve.  The  cyan  curve  denotes  the  residual  of  the 
terms  above  and  is  assumed  to  be  dominated  by  the  dissipation  of 
enstrophy  by  the  wavenumber  filter  in  the  model.  All  terms  non- 
dimensionalized  by  ( zq  -  u2)2/Ld. 


The  cyan  curve  again  denotes  the  imbalance  of  the  other 
terms.  Here,  the  imbalance  is  much  more  substantial  than 
in  the  energy  budget.  We  assume  that  the  imbalance  is 
dominated  by  small-scale  dissipation  (the  wavenumber 
filter),  which  is  much  more  important  in  enstrophy  bud¬ 
gets  than  in  energy  budgets.  As  in  Fig.  4,  we  see  in  Fig.  5 
that  the  fluxes  in  frequency  space  are  wider,  with  a  more 
rounded  appearance  than  they  have  in  wavenumber 
space. 

f  Necessity  for  averaging  over  many  grid  points 

It  is  necessary  to  average  over  many  grid  points  to 
obtain  an  estimate  of  the  spectral  flux  in  frequency  space. 
This  point  is  illustrated  in  Figs.  6  and  7.  Figure  6a  shows 
the  upper-layer  spectral  kinetic  energy  flux  in  wave- 
number  space  nKE,i(£),  computed  from  each  of  the  600 
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Fig.  6.  Variability  in  individual  realizations  of  spectral  fluxes  of 
upper-layer  kinetic  energy  nKE,i(k)  and  IIke,iG0  computed  from 
nominal  two-layer  FL  =  0.4  QG  simulation,  (a)  Black  curves  are 
I1ke,i(F)  computed  from  each  of  600  individual  snapshots  in  time, 
and  the  red  curve  is  nKE,i(k)  computed  from  averaging  over  the 
600  individual  snapshots,  (b)  Black  curves  are  nKE,i(w)  computed 
from  each  of  2562  individual  grid  points,  and  the  red  curve  is 
nKE,iH  computed  from  averaging  over  the  2562  grid  points.  All 
terms  nondimensionalized  by  (zq  -  u2)3 1  Ld: 


“realizations”  [individual  snapshots  (black  curves)]  of  the 
nominal  FL  =  0.4  QG  solution.  Averaging  over  the 
snapshots  (red  curve)  is  equivalent  to  integrating  over  all 
frequencies,  which  was  done  in  Fig.  4a.  Figure  6b  displays 
the  nKE,i(^)  plots  computed  at  the  2562  individual  grid 
points  in  the  QG  model  (black  curves),  alongside  the 
nKE,iH  values  computed  after  averaging  over  all  of 
these  grid  points  (red  curve;  equivalent  to  integrating 
over  all  wavenumbers,  as  in  Fig.  4b).  The  range  seen  in 
nKE,iH  estimates  at  individual  grid  points  is  about  25 
times  larger  than  the  range  seen  in  nKE,i(/c)  estimates  in 
individual  snapshots.  The  range  in  nKE,i(w)  values  is 
especially  large  for  low  frequencies. 

In  Figs.  7a, b,  we  show  the  spectral  fluxes  nKE,i(^)  and 
nKE,i(w),  respectively,  along  with  envelopes  represent¬ 
ing  the  variability  seen  in  Fig.  6.  The  envelopes  are 
computed  as  standard  deviations  of  the  individual  re¬ 
alizations  shown  in  Fig.  6,  divided  by  the  square  root  of 
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Fig.  7.  (a)  Upper-layer  spectral  kinetic  energy  flux  nKE,i(^) 
(solid  curve)  computed  from  nominal  two-layer  FL  =  0.4  QG 
simulation,  shown  along  with  envelopes  representing  variability  of 
individual  realizations  of  IIKE ^(k).  See  text  for  more  description, 
(b)  As  in  (a),  but  for  FIkeuC^)-  All  terms  nondimensionalized  by 
(Mi  ~  u2f/Ld. 


the  estimated  number  of  degrees  of  freedom.  The  num¬ 
ber  of  degrees  of  freedom  for  the  nKE,i(^)  calculation  is 
estimated  by  dividing  the  simulation  duration  by  the 
correlation  time  scale.  Similarly,  the  number  of  degrees 
of  freedom  for  the  nKE,i(^)  calculation  is  estimated  by 
dividing  the  number  of  grid  points  (25 62)  by  the  square  of 
the  correlation  length  scale.  Consistent  with  the  Fig.  6 
results,  the  nKE,i(^)  envelopes  are  much  smaller  than  the 
IIke,iH  envelopes. 

g.  Effects  of  simulation  duration  and  sampling 
frequency 

In  this  subsection,  we  examine  the  impact  of  the  du¬ 
ration  of  the  model  simulations  on  the  estimation  of 
spectra,  spectral  fluxes,  and  spectral  transfers  in  fre¬ 
quency  space.  We  are  motivated  in  part  by  the  large 
envelopes  computed  in  the  previous  subsection,  which 
suggest  a  large  uncertainty  in  the  estimation  of  nKE,i(^)* 
If  the  uncertainty  is  indeed  large,  then  different  simula¬ 
tions  of  the  nominal  FL  =  0.4  QG  run  should  yield  sub¬ 
stantially  different  estimates  of  nKE,i(<^)*  We  can  test  this 
with  simulations  having  different  durations.  We  will  ex¬ 
amine  spectral  transfers,  which  identify  sources  and  sinks 
of  energy.  We  are  interested  in  whether  spectral  trans¬ 
fers  will  identify  nonlinear  cascades  as  a  source  of  low- 
frequency  energy  and  whether  this  determination  will 
change  as  lower  frequencies  become  available  in  simu¬ 
lations  with  longer  durations. 

An  examination  of  a  nominal  FL  =  0.4  simulation  that 
is  saved  more  frequently  than  the  runs  discussed  thus  far 


will  shed  light  on  some  other  issues  of  interest.  We  an¬ 
ticipate  a  sensitivity  to  sampling  frequency  in  the  spectra 
at  high  frequencies  based  on  the  nature  of  small-scale 
dissipation.  Dissipation  operators  such  as  the  wave- 
number  filter  we  employ  here  act  preferentially  on  high 
wavenumbers  and  thus  yield  steep  slopes  at  the  high  end 
of  wavenumber  spectra  plots  (Figs.  2a, c).  The  wave- 
number  filter  does  not  act  directly  on  small  time  scales, 
and  therefore  it  is  not  clear  that  we  will  see  similarly 
steep  slopes  at  high  frequencies.  To  confidently  estimate 
the  slopes  at  high  frequencies,  we  need  simulations  that 
are  sampled  more  frequently.  Another  motivation  for 
examining  slopes  at  high  frequencies  was  noted  earlier. 
At  high  frequencies,  a  slope  steeper  than  -2  is  required 
for  the  first  moment  of  the  frequency  spectrum  coE  to  be 
a  meaningful  quantity. 

The  frequency  spectra  of  upper-layer  kinetic  energy  in 
three  simulations  of  the  nominal  FL  =  0.4  two-layer  QG 
solution,  with  each  having  a  different  duration  and  sam¬ 
pling  frequency,  are  displayed  in  Fig.  8a.  The  red  curve  is 
from  the  same  simulation  as  the  red  curves  in  Fig.  2:  600 
snapshots,  each  spaced  one  unit  of  nondimensional  time 
Ly(Zq  —  u2)  apart.  The  black  curve  is  also  sampled  one 
unit  of  nondimensional  time  apart  but  is  from  a  simula¬ 
tion  that  is  9  times  as  long  (5400  snapshots).  The  blue 
curve  is  from  a  simulation  that  contains  2400  snapshots, 
each  spaced  0.1  units  of  nondimensional  time  apart;  it  is 
more  frequently  sampled  but  shorter  than  the  other  two. 
The  space-  and  time -averaged  kinetic  energies  of  the 
three  simulations  are  equal  to  within  5%.  Over  the  win¬ 
dow  of  frequencies  utilized  in  Fig.  2b,  the  slope  of  the 
blue  curve  (higher-frequency  sampling)  is  —1.9,  similar  to 
the  slope  computed  in  that  band  for  the  red  curve  (—1.8). 
At  higher  frequencies,  the  slope  of  the  blue  curve  is 
steeper  than  -2,  as  anticipated  in  discussions  of  the  defi¬ 
nition  of  coE .  This  steeper  band  is  aliased  into  lower  fre¬ 
quencies  in  the  red  and  black  curves  (lower-frequency 
sampling)  in  Fig.  8a,  as  seen  by  the  upturns  at  the  high- 
frequency  end  of  those  curves. 

The  spectral  fluxes  of  upper-layer  kinetic  energy  in 
frequency  space  nKE1(cL>)  are  shown  in  Fig.  8b  for  the 
three  simulations.  Near  nondimensional  go  values  on  the 
order  of  one,  the  three  curves  agree  very  closely.  At  lower 
frequencies,  the  blue  curve  (short,  rapidly  sampled  sim¬ 
ulation)  separates  from  the  red  and  black  curves  by  slight 
amounts.  The  red  and  black  curves  (longer  simulations) 
stay  close  together  over  all  the  frequencies  for  which  both 
are  defined.  The  separation  of  the  three  curves  is  not 
nearly  as  large  as  the  envelopes  in  Fig.  7b  suggest  they 
might  be.  For  this  reason  we  prefer  the  term  “envelopes” 
in  the  discussion  of  Fig.  7b  rather  than  error  bars.  In  all 
three  curves  in  Fig.  8b,  the  nonlinear  kinetic  energy  flux 
extends  out  to  the  lowest  frequencies,  where  the  lowest 
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(a)  co  spectrum  of  (V  i^)2 


Fig.  8.  (a)  Frequency  spectra  and  (b)  spectral  fluxes  in  frequency 
space  nKE»  of  upper-layer  kinetic  energy  in  three  simulations 
of  the  nominal  two-layer  FL  =  0.4  QG  solution.  The  red  curve 
is  as  in  Fig.  2,  based  on  600  snapshots  spaced  one  unit  of  non- 
dimensional  time  Ldl(ux  ~u2)  apart.  The  blue  curve  is  based  on 
2400  snapshots,  spaced  0.1  units  of  nondimensional  time  apart.  The 
black  curve  is  based  on  5400  snapshots,  spaced  one  unit  of  non- 
dimensional  time  apart,  (c)  Spectral  transfers  T(a>)  =  -dU.(o))/do), 
shown  for  the  nominal  600-snapshot  simulation  (solid  curves)  and 
longer  5400-snapshot  simulation  (dashed  curves).  Colors  in  (c) 
represent  different  terms  in  the  model  energy  budget  as  in  Fig.  4. 


frequency  of  any  particular  run  is  set  by  the  duration  of 
the  run. 

As  discussed  in  many  sources  (e.g.,  Salmon  1978, 1980; 
Hua  and  Haidvogel  1986;  Larichev  and  Held  1995;  LaCasce 
1996;  among  others),  the  transfer  T(k )  =  -[dU(k)/dk] 
identifies  wavenumbers  that  serve  as  energy  sources  and 
sinks.  In  like  manner,  T(co)  =  -[dll(m)/dcu]  can  be  used 
to  identify  energy  sources  and  sinks  as  a  function  of 
frequency.  Figure  8c  shows  the  T(co)  budget  for  the 
600-snapshot  FL  =  0.4  QG  simulation  (solid  curves), 
next  to  the  budget  for  the  longer-duration  5400-snapshot 
simulation  (dashed  curves).  Different  terms  in  the  budget 
are  denoted  by  the  same  colors  as  in  Fig.  4.  At  low  fre¬ 
quencies,  the  main  balance  is  between  the  upper-layer 
kinetic  energy  flux  TKB1(co)  and  the  potential  energy  flux 


Fig.  9.  Spectral  flux  nKE,i(w)  °f  upper-layer  kinetic  energy  in 
frequency  space  computed  from  the  600-snapshot  nominal  two- 
layer  Fl  =  0.4  QG  simulation  (unfiltered)  and  from  filtered  ver¬ 
sions  of  this  simulation.  The  spatial  filter  is  a  Blackman  filter  that 
goes  to  zero  at  six  grid  points  from  the  central  grid  point.  The  time 
filter  is  a  Blackman  filter  that  goes  to  zero  at  four  units  of  non- 
dimensional  time  Ldl(u1  -  u2)  from  the  central  time.  The  wide  and 
very  wide  time  filters  are  Blackman  filters  that  go  to  zero  at  9  and  17 
units  of  nondimensional  time  from  the  central  time,  respectively. 
All  terms  nondimensionalized  by  (u1  —  u2)3 / 'Ld. 


7ape(&>),  with  substantial  contributions  from  forcing  and 
friction.  In  the  longer  run,  the  various  terms  in  the  budget 
extend  from  their  values  at  higher  frequencies  (which 
are  close  to  the  values  seen  in  the  shorter  run)  out  to 
lower  frequencies  in  a  relatively  flat  manner.  Figure  8c 
demonstrates  that  nonlinear  cascades  play  a  dominant 
role  in  maintaining  variance  out  to  frequencies  as  low 
as  one-thousandth  of  the  nondimensional  frequency 
(Uj  -  u2)/Ld.  We  have  not  yet  found  a  lower  bound  on 
the  frequencies  to  which  the  QG  turbulence  model  will 
push  into.  Considering  that  this  model  contains  only 
eddies,  which  are  normally  described  as  having  time 
scales  of  a  year  or  less,  and  has  no  gyres,  boundary 
currents,  or  other  features  that  are  normally  associated 
with  interannual  and  decadal  time  scales,  this  result  is 
somewhat  surprising. 

Given  that  the  kinetic  energy  flux  dominates  the  low- 
frequency  behavior,  it  is  of  interest  to  determine  how  large 
it  is  relative  to  the  total  energy  in  the  system.  The  absolute 
value  of  the  minimum  (trough)  of  the  flux  IIKE  i(m),  di¬ 
vided  by  the  domain-integrated  upper-layer  kinetic  energy 
in  the  simulation,  is  0.09,  in  units  of  (u{  —  u2)/Ld.  Taking 
typical  midocean  values  of  1-10  cm  s-1  for  ux  -  u2  and 
50  km  for  Ld. ,  this  amounts  to  2  X  10-8  to  2  X  10-7  s-1,  in 
dimensional  terms.  Later  we  will  compare  this  range  to 
values  computed  from  realistic  models  and  from  altimeter 
data. 
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Fig.  10.  Six  subdomains  used  to  compute  spectra  and  spectral  fluxes  in  V320  NLOM  output, 
highlighted  against  the  15  Feb  2002  snapshot  of  sea  surface  height  (cm)  in  the  model:  mid¬ 
latitude  southeast  Pacific,  high-latitude  southeast  Pacific,  Agulhas,  Malvinas,  Gulf  Stream,  and 
Kuroshio. 


h.  Impact  of  spatial  and  temporal  filtering  on  spectral 
fluxes  in  frequency  space 

In  Fig.  9,  motivated  by  the  spatial  and  temporal  fil¬ 
tering  inherent  in  the  creation  of  gridded  satellite  al¬ 
timeter  data,  we  show  nKE,i(^)  computed  from  the 
nominal  FL  =  0.4  QG  simulation  and  from  filtered  ver¬ 
sions  of  this  simulation.  Note  that  the  filtering  is  performed 
on  the  model  output  after  the  model  has  been  run,  not 
during  the  run.  The  time  filter  is  a  Blackman  filter  that  goes 
to  zero  at  four  units  of  nondimensional  time  Ldl(ju1  ~u2) 
from  the  central  time.  The  spatial  filter  goes  to  zero  at  six 
grid  points  from  the  central  grid  point,  meaning  it  has 
a  “full  width  at  half  maximum”  (width  at  the  half-power 
point)  of  3.4  grid  points  or  a  wavelength  of  0.83  Ld.  Fil¬ 
tering  in  space  or  in  time  diminishes  the  negative  lobe  of 
the  spectral  flux  substantially.  Filtering  in  both  space  and 
time  further  diminishes  the  negative  lobe.  Widening  the 
time  filter  distorts  the  spectral  flux  even  more  (cyan  curve), 
to  the  point  where  it  can  eventually  change  sign  (green 
curve).  We  will  revisit  this  point  later  in  the  paper. 

3.  Results  from  the  realistic  ocean  model 

We  use  three  full  years  of  output  (2001-03)  from  the 
nonassimilative  version  of  the  Naval  Research  Laboratory 


(NRL)  Layered  Ocean  Model  (NLOM).  Data-assimilative 
versions  of  NLOM  are  in  current  use  as  a  U.S.  Navy  op¬ 
erational  model.  NLOM  is  based  on  the  primitive  equa¬ 
tion  model  of  Hurlburt  and  Thompson  (1980),  but  with 
greatly  expanded  capability  (Wallcraft  et  al.  2003).  It  has 
six  dynamical  layers  and  a  bulk  mixed  layer.  Shriver  et  al. 
(2007)  and  references  therein  describe  the  NLOM  grid, 
wind  forcing,  and  other  model  details.  The  horizontal 
resolution  of  the  NLOM  simulation  analyzed  here  is  V320  in 
latitude  and  45/1024°  in  longitude,  on  a  model  grid  ex¬ 
tending  from  72°S  to  65°N.  For  simplicity,  we  refer  to  the 
simulations  as  having  V320  resolution.  The  simulation  we 
used  was  saved  at  6-h  intervals  over  the  3  yr  analyzed  here. 

Figure  10  displays  a  snapshot  of  sea  surface  height  in 
NLOM.  Six  subdomains  (regions)  of  interest,  in  which 
we  will  compute  spectra  and  spectral  fluxes,  are  shown. 
Four  of  the  subdomains  are  in  regions  of  intense  boundary 
current  and  mesoscale  eddy  activity:  the  Kuroshio  in  the 
western  North  Pacific,  the  Gulf  Stream  in  the  western 
North  Atlantic,  the  Agulhas  off  the  southern  tip  of  Af¬ 
rica,  and  the  Malvinas  in  the  western  South  Atlantic. 
Two  additional  subdomains  are  analyzed  in  the  eastern 
South  Pacific:  a  quiescent  region  in  midlatitudes  and 
a  high-latitude  region  that  touches  upon  the  Antarctic 
Circumpolar  Current  in  its  southernmost  portions. 
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In  the  subdomains  in  a  realistic  model  (or  observa¬ 
tional  oceanic  dataset),  spectral  flux  computations  will 
be  affected  by  energy  fluxes  into  and  out  of  the  sub- 
domains.  Scott  and  Wang  (2005)  investigated  the  impor¬ 
tance  of  subdomain  size  in  the  computation  of  spectral 
fluxes  nKE,i(£)  in  wavenumber  space.  They  computed 
fluxes  for  subdomains  10°,  20°,  and  40°  on  a  side.  The  10° 
subdomains  yielded  noisy  results,  “most  likely  indicating 
the  importance  of  boundary  terms  relative  to  interior 
terms.”  However,  their  Fig.  2  shows  that  even  the  noisy 
results  in  the  10°  X  10°  subdomains  are  not  qualitatively 
different  from  20°  X  20°  results  centered  on  the  same 
point,  which  in  turn  are  similar  to  the  40°  X  40°  results. 
Our  subdomains  in  this  paper  contain  roughly  as  many 
grid  points  as  the  20°  X  20°  boxes  in  their  paper.  There¬ 
fore,  our  spectral  fluxes  should  not  be  qualitatively  af¬ 
fected  by  the  size  of  our  subdomains. 

The  NLOM  computations  are  the  most  demanding 
calculations  shown  in  this  paper.  Three  years  of  6-hourly 
global  V320  NLOM  sea  surface  heights  represents  660  GB 
of  output  to  sort  through.  Once  extracted,  the  six  sub- 
domains  constitute  about  35  GB  of  output,  which  must  be 
analyzed  in  both  the  time  and  space  domains  several  times 
over  to  compute  nKE,i(<^)  and  nKE,i(^)*  For  simplicity, 
we  compute  only  IIke,i  in  the  NLOM  integral  energy 
budget  and  omit  other  terms.  Also  for  simplicity,  the 
factor  5/(1  +  8)  is  omitted  from  the  computation  of 
11ke,i>  because  stratification  is  not  uniform  across  the 
extent  of  the  model.  Finally,  the  envelopes  discussed  in 
section  2f  are  omitted  here  because  they  would  be  very 
expensive  to  compute. 

The  procedure  for  analyzing  NLOM  output  is  as  fol¬ 
lows:  The  NLOM  sea  surface  heights  are  first  low-pass 
filtered  in  time  to  remove  motions  with  periods  of  3  days 
or  less.  This  is  done  because  such  high-frequency  mo¬ 
tions  cannot  be  geostrophic.  Next,  we  compute  surface 
streamf unctions  t/q  =  giq/f,  where  17  is  sea  surface  height; 
g  =  9.8  m  s-2  is  the  gravitational  acceleration;  and /is 
the  Coriolis  parameter,  which  varies  as  the  sine  of  the 
model  gridpoint  latitude  (Vallis  2006).  Then,  for  each  of 
the  six  subdomains  shown  in  Fig.  10,  each  snapshot  of 
i/q  is  detrended  in  space  with  a  two-dimensional  least 
squares  fit.  Following  this,  the  i/q  field  is  tapered  in  space 
with  a  two-dimensional  function  constructed  of  nine 
overlapping  Hanning  windows  in  each  spatial  direction, 
with  each  window  having  a  width  one-fifth  of  the  domain 
extent  in  that  direction.  Tapering  is  also  done  on  each 
time  series  at  each  model  grid  point,  using  the  same  win¬ 
dow  as  in  the  QG  calculation  (nine  overlapping  Hanning 
windows).  This  type  of  detrending  and  tapering  is  standard 
in  computations  of  Fourier  transforms  of  nonperiodic 
data.  Fourier  transforms  are  then  repeatedly  utilized  to 
compute  derivatives.  For  simplicity,  we  base  the  increment 
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Fig.  11.  Frequency  spectra  of  the  variance  of  (a)  surface  ocean 
streamfunction  i/q  =  grj/f where  g  is  gravitational  acceleration, /is 
the  Coriolis  parameter,  and  17  is  sea  surface  height  drawn  from 
NLOM  output  in  subdomains  shown  in  Fig.  10,  and  (b)  geostrophic 
velocities  computed  from  this  streamfunction.  Dashed  lines  are 
drawn  at  the  frequencies  corresponding  to  periods  of  14  and  68.4 
days. 

lengths  of  the  wavenumbers  in  the  east-west  direction 
on  the  length  of  the  increments  at  the  central  latitude  of 
the  box.  Thus,  we  do  not  account  for  the  convergence  of 
meridians  in  our  spectral  calculations. 

The  frequency  spectra  of  the  variance  of  the  NLOM 
streamfunction  i/q  in  the  six  subdomains  are  displayed 
in  Fig.  11a.  Figure  lib  displays  the  frequency  spectra  of 
the  variance  of  geostrophic  velocities  (kinetic  energy; 
|Vi/q|2).  As  in  the  QG  results,  the  NLOM  frequency 
spectra  display  a  “flatness”  at  the  lowest  frequencies 
and  a  steeper  slope  at  higher  frequencies,  in  what  might 
be  called  the  “mesoscale  band.”3  The  dashed  vertical 


3  Schmitz  (1996)  (see  also  references  therein,  including  Schmitz 
and  Luyten  1991)  defines  the  mesoscale  band  as  20-200  days.  This 
band  contains  most  of  the  kinetic  energy  in  regions  of  strong 
boundary  currents. 
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Fig.  12.  Spectral  fluxes  of  surface  ocean  geostrophic  kinetic  energy  IIke,i(0  in  wavenumber 
space  computed  from  NLOM  output  in  subdomains  shown  in  Fig.  10. 


lines  in  Fig.  11  are  at  frequencies  corresponding  to  pe¬ 
riods  of  68.4  and  14  days,  respectively.  In  the  next  sec¬ 
tion,  we  will  discuss  the  slopes  in  this  band  along  with 
slopes  computed  in  the  same  band  from  satellite  altimeter 
data.  This  band  was  chosen  because  it  lies  in  the  meso- 
scale  band,  is  resolved  by  both  the  NLOM  and  AVISO 
results  utilized  in  this  paper,  appears  to  lie  in  approxi¬ 
mately  the  same  part  of  parameter  space  in  which  the 
band  analyzed  in  Figs.  2b, d  lies  (i.e.,  in  the  steep  part  of 
the  QG  spectrum),  and  covers  a  range  of  frequencies 
equal  to  the  range  of  the  band  analyzed  in  Figs.  2b, d.  The 
highest  frequencies  in  the  NLOM  spectra  show  a  steep 
dropoff  due  to  the  3 -day  low-pass  filter  we  employed. 

Figure  12  displays  the  spectral  fluxes  nKE,i(^)  in 
wavenumber  space  computed  from  the  six  subdomains 
in  NLOM.  As  in  previous  studies  (Scott  and  Wang  2005; 


Scott  and  Arbic  2007;  Schlosser  and  Eden  2007;  Tulloch 
et  al.  2011),  the  negative  lobes  representing  inverse  cas¬ 
cades  are  ubiquitous  (seen  in  every  subplot).  As  noted  in 
Scott  and  Wang  (2005),  the  magnitude  of  the  spectral 
fluxes  varies  by  two  orders  of  magnitude  across  ocean 
basins. 

Figure  13  displays  the  spectral  fluxes  nKE,i(w)  in  fre_ 
quency  space  (shown  as  thick  solid  curves)  computed 
from  the  six  subdomains  in  NLOM.  The  fluxes  are  dom¬ 
inated  by  negative  lobes,  implying  that  nonlinear  terms 
drive  a  cascade  of  kinetic  energy  toward  longer  time 
scales.  Comparison  of  the  subplots  in  Fig.  13  to  the  sub¬ 
plots  in  Fig.  12  demonstrates  that  the  magnitudes  of  the 
spectral  fluxes  nKE,i(w)  in  the  various  subdomains  are 
comparable  to  but  smaller  than  the  magnitudes  of  the 
corresponding  spectral  fluxes  nKE,i(^)*  The  ratio  of  the 
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(a)  Mid-latitude  SE  Pacific  (b)  High-latitude  SE  Pacific 


Fig.  13.  Thick  solid  curves  are  spectral  fluxes  of  surface  ocean  geostrophic  kinetic  energy 
nKE,iW  in  frequency  space  computed  from  NLOM  output  in  subdomains  shown  in  Fig.  10. 
Thick  dashed  curves  represent  fluxes  computed  after  the  NLOM  output  has  been  temporally 
filtered.  See  text  for  details. 


minimum  (trough)  of  the  flux  I1KE to  the  domain- 
averaged  kinetic  energy  lies  in  the  range  from  5.8  X  10-8 
to  1.4  X  10-7  s^1,  similar  to  the  range  estimated  from 
the  QG  results  in  section  2g.  This  suggests  to  us  that,  as 
in  the  QG  simulations,  nonlinearities  play  an  important 
role  in  the  maintenance  of  low-frequency  energy  in  NLOM. 

The  dashed  thick  curves  in  Fig.  13  are  spectral  fluxes 
nKE,iW  computed  from  NLOM  output  that  has  been 
temporally  smoothed,  in  an  attempt  to  mimic  the  fil¬ 
tering  inherent  in  the  gridded  satellite  altimeter  data.  As 
in  the  calculations  made  from  filtered  QG  model  output, 
the  filtering  is  performed  on  NLOM  output  after  the 
model  has  been  run,  not  during  the  runs.  For  NLOM  we 
use  a  Blackman  filter,  which  goes  to  zero  at  7  days  from 


the  central  time.  The  7-day  width  is  chosen  to  match  the 
7-day  output  frequency  of  the  AVISO  data.  As  in  the 
QG  results,  temporal  smoothing  diminishes  the  inverse 
temporal  cascade,  to  the  point  where  in  some  cases  (see 
especially  Fig.  13e)  the  flux  changes  sign. 

4.  Results  from  satellite  altimeter  data 

We  use  840  snapshots  of  the  two-satellite  AVISO  V30 
Mercator  grid  “reference”  product  (Le  Traon  et  al. 
1998;  Ducet  et  al.  2000),  beginning  with  14  October  1992 
and  ending  with  12  November  2008.  The  AVISO  record 
we  utilize  is  much  longer  than  the  NLOM  record  (16 
versus  3  yr),  but  has  a  much  lower  temporal  sampling 
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Fig.  14.  Six  subdomains  used  to  compute  spectra  and  spectral  fluxes  in  lh°  AVISO  gridded 
satellite  altimeter  data,  highlighted  against  the  27  Dec  2006  snapshot  of  sea  surface  height 
anomaly  (cm)  in  AVISO:  midlatitude  southeast  Pacific,  high-latitude  southeast  Pacific,  Agulhas, 
Malvinas,  Gulf  Stream,  and  Kuroshio. 


rate  (7  days  versus  6  h).  The  spatial  resolution  is  also 
much  lower  (V30  versus  V320).  The  analysis  procedure 
used  on  the  AVISO  data  is  very  similar  to  that  used  on 
the  NLOM  output,  with  the  main  exception  being  that 
we  do  not  subject  the  AVISO  data  to  the  3-day  low-pass 
filtering. 

Figure  14  displays  the  six  subdomains  over  which  we 
compute  spectra  and  spectral  fluxes  from  AVISO  data. 
The  subdomains  correspond  closely  to  those  utilized  in 
the  NLOM  analysis.  The  correspondence  is  not  exact 
because  of  the  lack  of  points  in  shallow  waters  in  NLOM 
(see  Shriver  et  al.  2007,  and  references  therein).  The 
misfit  between  the  boundaries  of  the  NLOM  and  altim¬ 
eter  subdomains  is  largest  in  the  Agulhas  subdomain.  As 
in  the  NLOM  computation,  we  compute  only  nKE,i  from 
AVISO,  omitting  other  terms  in  the  integral  energy 
budget,  and  we  omit  the  stratification  factor  5/(1  +  5). 

Figure  15  displays  the  frequency  spectra  of  the  vari¬ 
ance  of  i/q  and  of  surface  geostrophic  velocity  computed 
from  AVISO.  As  in  the  NLOM  and  idealized  QG  results, 
the  AVISO  spectra  are  flat  at  the  lowest  frequencies  and 
show  a  steeper  slope  at  higher  frequencies. 

Table  1  presents  the  slopes  of  the  frequency  spectra  in 
NLOM  and  AVISO,  over  the  period  band  14-68.4  days, 
for  the  six  subdomains  and  for  both  streamfunction  and 


velocity  variance.  A  few  general  conclusions  can  be 
noted.  For  the  most  part,  the  NLOM  spectral  slopes 
hover  around  the  —1.8  and  -2.6  values  found  in  the 
high-frequency  QG  model  results,  albeit  with  wide 
variability  (deviations  vary  by  as  much  as  0.4  in  either 
direction).  A  notable  exception  is  the  slope  of  the 
streamfunction  variance  in  the  midlatitude  southeast 
Pacific,  which  is  actually  less  steep  than  the  slope  of  the 
velocity  variance  (at  higher  frequencies,  the  stream- 
function  variance  in  this  subdomain  is  steeper  than  the 
velocity  variance,  as  is  more  typically  the  case).  Con¬ 
sistent  with  the  highly  smoothed  nature  of  the  AVISO 
product,  the  AVISO  slopes  are  steeper  than  the  NLOM 
slopes,  often  substantially  so  (deviations  vary  by  as  much 
as  2.0).  With  the  exception  of  the  midlatitude  southeast 
Pacific  patch  noted  above,  the  difference  between  slopes 
in  NLOM  streamfunction  variance  and  NLOM  velocity 
variance  spectra,  or  in  AVISO  streamfunction  variance 
and  AVISO  velocity  variance  spectra,  ranges  from  -0.3 
to  -0.7.  These  differences  are  remarkably  consistent  be¬ 
tween  AVISO  and  NLOM  output  computed  for  the  same 
subdomain,  despite  the  fact  that  the  AVISO  slopes  are 
generally  much  steeper  than  the  NLOM  slopes. 

The  values  of  the  slopes  depend,  of  course,  on  the  band 
chosen  for  computing  the  least  squares  fit.  The  slopes  of 
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(a)  co  spectrum  of  \| 


Fig.  15.  Frequency  spectra  of  the  variance  of  (a)  surface  ocean 
streamfunction  i f/t  =  grjlf,  where  g  is  gravitational  acceleration, /is 
the  Coriolis  parameter,  and  17  is  sea  surface  height  drawn  from 
gridded  V30  AVISO  data  in  subdomains  shown  in  Fig.  14,  and  (b) 
geostrophic  velocities  computed  from  this  streamfunction.  Dashed 
lines  are  drawn  at  the  frequencies  corresponding  to  periods  of  14 
and  68.4  days. 

the  AVISO  spectra  are  not  very  clean  in  the  14-68.4-day 
band.  Note,  for  instance,  the  upturns  at  the  high-frequency 
end  of  the  Kuroshio  results  in  Fig.  15.  This  is  likely  due 
to  aliasing;  see  the  discussion  of  the  upturns  in  Fig.  8a. 
We  also  computed  spectral  slopes  over  the  38-184-day 


band,  which  gets  us  away  from  the  problems  at  the  high- 
frequency  end  for  AVISO.  For  the  sake  of  brevity,  the 
slopes  over  the  38-184-day  band  are  not  given  here  in 
detail,  but  we  do  note  that  over  most  of  the  six  sub- 
domains  they  are  more  consistent  between  NLOM  and 
AVISO  than  the  14-68.4-day  band  results.  This  does  not 
mean  that  the  38-184-day  band  is  without  problems.  It  is 
not  as  cleanly  separated  from  the  flat  low-frequency  part 
of  the  spectrum  as  the  14-68.4-day  band  is.  Thus,  the 
slopes  in  the  38-184-day  band  are  more  likely  to  rep¬ 
resent  a  mixture  of  the  dynamics  of  two  distinct  regimes. 
The  38-184-day  band  slopes  are  generally  shallower 
than  the  14-68.4-day  band  slopes. 

Periods  shorter  than  14  days  are  not  represented  in  the 
AVISO  data,  whereas  periods  longer  than  1095  days  are 
not  represented  in  the  3-yr  NLOM  record  used  here. 
Slopes  computed  from  NLOM  output  over  a  frequency 
band  corresponding  to  periods  between  3  and  14.6  days 
are  considerably  steeper  than  those  in  Table  1.  For  the 
3-14.6-day  band,  the  differences  in  slopes  between  NLOM 
streamfunction  and  velocity  variance  spectra  range  from 
-0.5  to  -1.1,  except  in  the  high-latitude  southeast  Pa¬ 
cific  patch,  which  yields  a  slope  that  is  steeper  in  velocity 
variance  by  0.2.  For  periods  between  294  and  2940  days, 
the  slopes  of  the  AVISO  frequency  spectra  range  from 
-0.5  to  -0.1  for  streamfunction  variance  and  from  -0.2 
to  0.1  for  velocity  variance,  with  the  streamfunction 
variance  always  being  steeper,  by  amounts  ranging  from 
-0.02  to  -0.3.  These  values  lie  near  the  low-frequency 
QG  values  reported  earlier  of  -0.4  and  -0.1  for  stream- 
function  and  velocity  variance  and  -0.2  difference  be¬ 
tween  these. 

All  six  of  the  AVISO  subdomains  display  a  tendency 
for  the  low-frequency  variance  of  zonal  velocity  to  be 
larger  than  that  of  meridional  velocity.  Figure  16  plots 
the  frequency  spectra  of  zonal  and  meridional  velocity 
variance  in  the  midlatitude  southeast  Pacific,  in  which 
this  tendency  is  most  extreme.  The  zonal  velocity  vari¬ 
ance  integrated  over  periods  of  420-5880  days  is  2.3 
times  larger  than  the  meridional  velocity  variance.  In  the 


Table  1.  The  first  column  is  the  six  subdomains  of  World  Ocean  examined  in  both  NLOM  output  and  AVISO  gridded  satellite 
altimeter  data.  The  second  through  fifth  columns  are  the  slopes  of  frequency  spectra  of  variance  of  surface  streamfunction  i/q  and  variance 
of  surface  geostrophic  velocity  (kinetic  energy;  |  Vi/q|2).  Slopes  are  computed  as  least  squares  fits  over  the  frequency  band  corresponding  to 
periods  of  14-68.4  days  (see  dashed  vertical  lines  in  Figs.  11, 15).  The  sixth  column  is  NLOM  A,  the  slope  in  the  NLOM  streamfunction 
variance  spectrum  minus  the  slope  in  the  NLOM  velocity  variance  spectrum.  The  seventh  column  is  as  in  the  sixth  column,  but  for  AVISO. 
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Fig.  16.  Frequency  spectra  of  the  variance  of  zonal  and  meridi¬ 
onal  geostrophic  velocity,  computed  from  AVISO  gridded  altimeter 
data  in  the  midlatitude  southeast  Pacific  patch  shown  in  Fig.  14. 


other  five  subdomains,  the  zonal  velocity  variance  in¬ 
tegrated  over  these  periods  ranges  from  1.1  to  1.6  times 
larger  than  the  meridional  variance.  This  tendency  for 
zonal  velocity  variance  to  exceed  meridional  velocity 
variance  at  low  frequencies  was  noted  by  Richman  et  al. 
(1977),  from  analysis  of  a  current-meter  record  much 
shorter  than  the  present-day  AVISO  record. 

Figure  17  displays  the  spectral  fluxes  nKE,i(^)  in 
wavenumber  space  over  the  six  subdomains  of  AVISO 
data.  As  in  the  QG  and  NLOM  computations,  the  dom¬ 
inant  feature  is  a  negative  lobe,  representing  an  inverse 
cascade  toward  larger  length  scales.  Also  as  in  the  NLOM 
computations,  the  magnitudes  of  the  spectral  fluxes  vary 
by  two  orders  of  magnitude  over  ocean  basins.  Envelopes 
were  computed  assuming  a  35-day  correlation  time  over 
the  16-yr  record,  and  are  small  enough  that  they  are  dif¬ 
ficult  to  see  in  Fig.  17. 

Figure  18  displays  the  spectral  fluxes  nKE,i(w)  in  fre_ 
quency  space  over  the  six  subdomains  of  AVISO  data. 
To  test  the  sensitivity  to  windowing,  fluxes  are  computed 
in  three  different  ways.  The  baseline  calculations  (black 
curves)  use  the  full  time  interval  of  840  snapshots,  ta¬ 
pered  with  nine  overlapping  Hanning  windows  each  168 
snapshots  wide.  We  also  subdivided  the  time  series  into 
time  intervals  of  160  days  each  with  9  different  over¬ 
lapping  realizations  (blue  curves)  and  into  time  intervals 
of  80  days  each  with  19  different  overlapping  realizations 
(red  curves).  At  higher  frequencies,  the  three  methods 
yield  similar  results,  suggesting  that  our  computation 
technique  is  relatively  robust.  Some  subdomains  (high- 
latitude  southeast  Pacific,  Agulhas,  and  arguably  the 
Malvinas)  display  a  negative  lobe  indicating  an  inverse 


temporal  cascade,  consistent  with  the  QG  and  NLOM 
results.  However,  the  Gulf  Stream,  Kuroshio,  and  mid¬ 
latitude  southeast  Pacific  display  a  predominantly  forward 
cascade  (positive  spectral  flux)  in  frequency  space.  This 
may  be  due  to  interesting  physics  in  some  regions  of  the 
ocean  that  are  not  captured  well  in  either  NLOM  or  the 
QG  turbulence  model.  The  smoothing  exercises  in  pre¬ 
vious  sections  suggest  an  alternative  possibility,  that  for¬ 
ward  temporal  cascades  could  be  artifacts  of  the  filtering 
inherent  in  the  creation  of  gridded  satellite  altimeter  data 
products.  For  the  subdomains  with  a  clear  negative  spec¬ 
tral  flux  (Agulhas  and  high-latitude  southeast  Pacific),  the 
ratios  of  the  minimum  (trough)  negative  flux  nKE,i(w) to 
the  kinetic  energy  are  6.0  and  4.0  X  10-8  s-1,  respectively, 
consistent  with  the  QG  and  NLOM  results  (see  discussions 
in  earlier  sections)  and  suggestive  of  the  importance  of 
nonlinearities  in  maintaining  energy  at  low  frequencies. 

5.  Summary  and  discussion 

Motivated  by  the  ubiquity  of  time  series  in  oceanic 
data,  the  relative  lack  of  studies  of  geostrophic  turbu¬ 
lence  in  the  frequency  domain,  and  the  interest  in  quan¬ 
tifying  the  contributions  of  intrinsic  oceanic  nonlinearities 
to  oceanic  frequency  spectra,  we  have  undertaken  a  com¬ 
prehensive  examination  of  spectra  and  spectral  fluxes  of 
surface  oceanic  geostrophic  flows  in  the  frequency  do¬ 
main.  We  have  shown  that  the  frequency  spectra  of  the 
streamfunction  (sea  surface  height)  variance  and  geo¬ 
strophic  velocity  variance  (kinetic  energy)  in  idealized 
flat-bottom  f-plane  two-layer  QG  turbulence  models, 
realistic  ocean  general  circulation  models,  and  gridded 
satellite  altimeter  data  exhibit  qualitative  similarities, 
with  steeper  slopes  at  higher  frequencies  and  flatter  slopes 
at  lower  frequencies.  In  the  data  and  models  considered 
here,  the  flat  slopes  at  lower  frequencies  continue  out 
to  the  lowest  frequencies  allowed  by  the  duration  of  the 
record. 

We  demonstrate  here  that  the  machinery  of  geo¬ 
strophic  turbulence,  which  is  usually  discussed  in  the 
wavenumber  domain,  is  useful  in  the  frequency  domain 
as  well.  Through  examination  of  the  energy  budgets  of 
the  QG  model  in  frequency  space,  we  find  that  the 
quantity  nKE,i(w)>  the  spectral  flux  of  upper-layer  ki¬ 
netic  energy  in  the  frequency  domain,  documents  an 
energy  cascade  toward  longer  time  scales  driven  by  the 
action  of  nonlinearities,  much  like  nKE,i(^)>  the  spectral 
flux  of  upper-layer  kinetic  energy  in  the  wavenumber 
domain,  documents  an  inverse  cascade  of  energy  toward 
larger  length  scales.  The  IIke,^)  is  “wider”  than 
11ke,i(^)>  consistent  with  the  action  of  nonlinearities  at 
the  longest  time  scales  permitted  by  the  record  duration. 
A  budget  of  the  spectral  transfers  T(co)  =  —dll (co)/dco 
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(a)  Midlatitude  SE  Pacific 


(c)  Agulhas 


(e)  Gulf  Stream 
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(b)  High  latitude  SE  Pacific 
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Fig.  17.  Spectral  fluxes  nKE,i(0  in  wavenumber  space  (thick  curves)  computed  from  gridded 
V30  AVISO  data  in  subdomains  shown  in  Fig.  14.  Envelopes  are  drawn  with  dashed  curves. 


further  demonstrates  the  importance  of  nonlinearities  in 
maintaining  low-frequency  energy  in  the  QG  turbulence 
model  examined  here. 

We  have  also  calculated  spectral  fluxes  in  frequency 
space  from  NLOM,  a  realistic  near-global  ocean  general 
circulation  model,  and  from  the  AVISO  gridded  satellite 
altimeter  dataset.  The  IIke,^)  values  computed  from 
the  realistic  model  also  document  a  nonlinear  cascade  of 
kinetic  energy  toward  longer  time  scales.  The  IIke,:^) 
values  computed  from  AVISO  gridded  satellite  altimeter 
data  yield  more  ambiguous  results.  Some  subdomains 
(regions)  are  dominated  by  an  “inverse  temporal  cas¬ 
cade,”  whereas  others  exhibit  a  forward  temporal  cas¬ 
cade.  Computations  done  with  temporally  and/or  spatially 
smoothed  QG  model  output,  as  well  as  temporally 


smoothed  NLOM  output,  demonstrate  that  the  spectral 
flux  computations  are  highly  susceptible  to  the  smoothing 
inherent  in  the  construction  of  gridded  altimeter  prod¬ 
ucts.  It  remains  to  be  seen  whether  the  forward  temporal 
cascades  seen  in  some  subdomains  in  altimeter  data 
represent  real  physics  not  present  in  the  models  analyzed 
here  or  artifacts  of  the  smoothed  nature  of  gridded  al¬ 
timeter  data. 

A  comprehensive  comparison  of  the  slopes  of  fre¬ 
quency  spectra  in  the  QG  turbulence  model,  the  realistic 
NLOM  simulation,  and  AVISO  altimetry  data  is  given 
in  this  paper.  In  the  QG  model,  the  slope  of  the  velocity 
variance  spectra  at  the  high-frequency  end  is  —1.8,  close 
to  the  —2  prediction  from  dimensional  considerations  of 
a  kinetic  energy  cascade  in  frequency  space  (Tennekes 
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(a)  Midlatitude  SE  Pacific  (b)  High  latitude  SE  Pacific 
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Fig.  18.  Spectral  fluxes  nKE,i(w)  in  frequency  space  (thick  curves)  computed  from  gridded 
V30  AVISO  data  in  subdomains  shown  in  Fig.  14.  Thick  black  curves  denote  baseline  results 
computed  over  840  weekly  snapshots.  Thick  blue  (red)  curves  denote  results  computed  over 
overlapping  intervals  of  160  (80)  weekly  snapshots  each.  Dashed  curves  denote  envelopes  of 
thick  black  curves. 


and  Lumley  1972).  The  slopes  in  frequency  spectra  of 
streamf unction  and  velocity  variance  do  not  differ  by 
-2,  as  they  do  in  wavenumber  spectra,  which  represents 
a  quantitative  failure  of  the  Taylor  (1938)  hypothesis 
relating  frequency  to  wavenumber  spectra.  In  the  band 
of  higher  frequencies  analyzed  in  the  QG  model,  the 
slopes  differ  by  -0.8.  In  the  band  14-68.4  days,  the  dif¬ 
ferences  in  the  slopes  of  the  streamfunction  and  velocity 
variance  spectra  vary  between  -0.3  and  —0.7  in  NLOM 
and  AVISO,  in  five  of  the  six  subdomains  analyzed. 


There  is  some  consistency  between  slopes  in  the  ideal¬ 
ized  QG  model  and  the  realistic  model  (NLOM),  but  the 
latter  shows  wide  variability  from  one  region  to  another. 
In  the  14-68.4-day  band,  the  slopes  in  the  AVISO  dataset 
are  significantly  steeper,  consistent  with  the  highly  filtered 
nature  of  that  product.  The  failure  of  the  Taylor  (1938) 
hypothesis  in  our  results  is  reminiscent  of  discussions  in 
Hayashi  and  Golder  (1977),  who  found  no  simple  re¬ 
lationship  between  kinetic  energy  spectra  in  the  fre¬ 
quency  and  wavenumber  domains  in  their  space-time 
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analysis  of  realistic  atmospheric  models  and  atmospheric 
datasets. 

The  time  scales  in  two-layer  QG  turbulence  exhibit 
a  sensitivity  to  the  strength  of  linear  bottom  Ekman 
friction,  somewhat  similar  to  the  sensitivity  previously 
documented  in  the  length  scales.  Although  there  are 
some  interesting  differences  discussed  in  section  2c,  an 
important  similarity  is  that  the  time  scales  in  two-layer 
QG  turbulence  tend  to  be  longer  than  the  imposed  non- 
dimensional  time  Ldl(jT[  -  u ^),  just  as  the  length  scales 
tend  to  be  larger  than  the  imposed  length  scale  Ld. 

We  end  by  briefly  discussing  natural  extensions  of  this 
work.  In  this  study,  we  have  looked  only  at  the  spectra 
and  spectral  fluxes  of  geostrophic  flows.  In  future  work, 
we  will  consider  spectra  and  spectral  fluxes  computed 
from  the  full  velocity  field  in  realistic  models  that  con¬ 
tain  ageostrophic  flows.  In  this  study,  we  have  looked 
only  at  spectra  and  spectral  fluxes  in  wavenumber  space 
or  in  frequency  space,  but  never  in  mixed  frequency- 
wavenumber  space.  In  section  2d,  we  developed  the 
theory  for  the  more  general  flux  H(k,  co).  We  anticipate 
that  examination  of  spectra  in  mixed  wavenumber- 
frequency  space  may  help  to  illuminate  some  of  the 
limitations  of  Taylor’s  (1938)  hypothesis.  We  antici¬ 
pate  that  examination  of  mixed  fluxes  II (/c,  co),  in  flows 
that  include  ageostrophic  motions  such  as  near-in- 
ertial  waves  and  tides  as  well  as  geostrophic  eddies, 
could  yield  insights  into  interactions  between  these 
different  classes  of  motions.  Motivated  by  the  results 
in  Fig.  8c,  we  also  anticipate  that  the  spectral  transfers 
T(co )  could  be  useful  diagnostics  of  the  role  of  non¬ 
linearity  in  maintaining  low-frequency  variability  in  the 
climate  system.  In  ocean  models,  T(co)  could  be  com¬ 
puted  for  both  the  atmospheric  forcing  and  the  oceanic 
nonlinearities,  potentially  helping  to  determine  how 
much  low-frequency  variability  in  the  ocean  is  driven  by 
atmospheric  forcing  (e.g.,  Griffies  and  Tziperman  1995) 
and  how  much  is  due  to  intrinsic  oceanic  nonlinearities 
(e.g.,  Penduff  et  al.  2011). 
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